/* compute pi using Monte Carlo method */
#include "mpi.h"
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define CHUNKSIZE 1000
/* message tags */
#define REQUEST 1
#define REPLY 2

int main(int argc, char *argv[]) {
  int iter;
  int in, out, i, iters, max, ix, iy, ranks[1], done, temp;
  double x, y, Pi, error, epsilon;
  int numprocs, myid, server, totalin, totalout, workerid;
  int rands[CHUNKSIZE], request;
  MPI_Comm world, workers;
  MPI_Group world_group, worker_group;
  MPI_Status status;

  MPI_Init(&argc, &argv);
  world = MPI_COMM_WORLD;
  MPI_Comm_size(world, &numprocs);
  MPI_Comm_rank(world, &myid);
  server = numprocs - 1; /* last proc is server */
  if (myid == 0) {
    if (argc < 2) {
      fprintf(stderr, "Usage: %s epsilon\n", argv[0]);
      MPI_Abort(MPI_COMM_WORLD, 1);
    }
    sscanf(argv[1], "%lf", &epsilon);
  }
  MPI_Bcast(&epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Comm_group(world, &world_group);
  ranks[0] = server;
  MPI_Group_excl(world_group, 1, ranks, &worker_group);
  MPI_Comm_create(world, worker_group, &workers);
  MPI_Group_free(&worker_group);
  if (myid == server) { /* I am the rand server */
    do {
      MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST, world, &status);
      if (request) {
        for (i = 0; i < CHUNKSIZE;) {
          rands[i] = random();
          if (rands[i] <= INT_MAX)
            i++;
        }
        MPI_Send(rands, CHUNKSIZE, MPI_INT, status.MPI_SOURCE, REPLY, world);
      }
    } while (request > 0);
  } else { /* I am a worker process */
    request = 1;
    done = in = out = 0;
    max = INT_MAX; /* max int, for normalization */
    MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
    MPI_Comm_rank(workers, &workerid);
    iter = 0;
    while (!done) {
      iter++;
      request = 1;
      MPI_Recv(rands, CHUNKSIZE, MPI_INT, server, REPLY, world,
               MPI_STATUS_IGNORE);
      for (i = 0; i < CHUNKSIZE;) {
        x = (((double)rands[i++]) / max) * 2 - 1;
        y = (((double)rands[i++]) / max) * 2 - 1;
        if (x * x + y * y < 1.0)
          in++;
        else
          out++;
      }
      MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM, workers);
      MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM, workers);
      Pi = (4.0 * totalin) / (totalin + totalout);
      error = fabs(Pi - 3.141592653589793238462643);
      done = (error < epsilon || (totalin + totalout) > 100000000);
      request = (done) ? 0 : 1;
      if (myid == 0) {
        printf("\rpi = %23.20f", Pi);
        MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
      } else {
        if (request)
          MPI_Send(&request, 1, MPI_INT, server, REQUEST, world);
      }
    }
    MPI_Comm_free(&workers);
  }

  if (myid == 0) {
    printf("\npoints: %d\nin: %d, out: %d, <ret> to exit\n", totalin + totalout,
           totalin, totalout);
    getchar();
  }
  MPI_Finalize();
  return 0;
}